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Abstract. Force-based atomistic-continuum hybrid methods are the only known pointwise 
consistent methods for coupling a general atomistic model to a finite element continuum 
model. For this reason, and due to their algorithmic simplicity, force-based coupling methods 
have become a popular class of atomistic-continuum hybrid models as well as other types 
of multiphysics models. However, the recently discovered unusual stability properties of the 
linearized force-based quasicontinuum (QCF) approximation, especially its indefiniteness, 
present a challenge to the development of efficient and reliable iterative methods. 

We present analytic and computational results for the generalized minimal residual (GM- 
RES) solution of the linearized QCF equilibrium equations. We show that the GMRES 
method accurately reproduces the stability of the force-based approximation and conclude 
that an appropriately preconditioned GMRES method results in a reliable and efficient 
solution method. 



1. Introduction 

The motivation for coupled atomistic/continuum models of solids is that the accuracy of 
an atomistic model is often only needed in localized regions of the computational domain, but 
a coarse-grained continuum model is necessary for the simulation of large enough systems to 
include long-ranee effects The force-based approach has become 

very popular because it provides a particularly simple and accurate (T3J method for coupling 
two physics models without the development of an accurate hybrid coupling energy. It 
operates by creating disjoint subdomains in which the equilibrium equations at each degree 
of freedom are obtained by assigning forces directly from one of the physics models. In 
addition to coupling atomistic and continuum models, such an approach has also been found 
to be attractive, for example, in the coupling of regions modeled by quantum mechanics 
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to regions modeled by molecular mechanics, since accurate hybrid coupling energies require 
an interfacial region that is too computationally demanding for the quantum mechanics 
model jl]. 

The force-based quasicontinuum (QCF) approximation is attractive because of its simple 
and efficient implementation and because it is the only known pointwise consistent quasi- 
continuum (QC) approximation for coupling a general atomistic model with a Cauchy-Born 
continuum model [13]. By consistent we mean that the absence of ghost forces under ho- 
mogeneous deformations. Its main drawback is that it results in a non-conservative force 
field j6], that is, the QCF forces are not compatible with any energy functional. Several cre- 
ative attempts have been made to develop hybrid coupling energies that satisfy the patch test 



(there are no resultant forces under uniform strain) [14 29 , which is a weaker compatibility 
condition than pointwise consistency and leads to reduced accuracy. 

In this paper, we consider the force-based quasicontinuum approximation (QCF), 

-^ qc V cf ) = z, (i) 

but, for simplicity, we will focus mainly on its linearization about a reference state, 

Lfu^ = /; 

sec Section [2] for the precise definitions. Recent analyses of the linearized QCF opera- 
tor |l2"p3"] have identified both further advantages as well as disadvantages of the force-based 
coupling approach. In addition to being non-symmetric, which is related to the fact that 
jrqct j g non _ conservative, the linearized QCF operator also suffers from a lack of positive- 



definiteness 13 . In the present paper, we show that this somewhat unusual stability property 



of the operator L^ f presents a challenge for the development of efficient and stable iterative 
solution methods that is overcome by the GMRES methods we propose. 

1.1. Framework for iterative solution methods. We consider three related approaches 
to the development of iterative methods for the QCF equilibrium equations Q. A pop- 
ular approach [2l] to solve the force-based equations Q modifies a nonlinear conjugate 
gradient algorithm by replacing the univariate optimization of an energy, used for step size 



selection 23 , with the computation of a step size such that the residual is (approximately) 
orthogonal to the current search direction. We will show in Section |4.2| that, due to the 
indefiniteness of L^ cf , this method is not numerically stable for our QCF model problem. 
The second approach we consider is the nonlinear splitting 

-^ c \y) = - [^ cf (y) + V£(y)] + VS(y) 

to construct the nonlinear iteration equation 

V£(yW) = f+ [j- qcf (yW) + V£(y<">)] . (2) 

The iterative solution of the nonlinear splitting method ^ can then be obtained from the 
minimization of the sum of S(y) and the potential energy of the dead load / + where 

gin) . = J n* („(»)) +V £ ( y W) j 

that is, 

G argmin {y ^ E{y) - (f,y) - (g^,y)}. 

For this approach to be accurate under conditions near the formation or motion of defects, 
care must be taken to ensure that the energy E(y) accurately reproduces the stability of the 
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approximated atomistic system. We will see in Section 4.1 that using the original quasi 



continuum energy £ qce (y) defined in (18), which results in the ghost force correction (GFC) 



scheme, does not reliably reproduce the stability of the atomistic system 10 and can give a 
reduced critical strain for a lattice instability. 

To develop the final approach, we recall the Newton method 

- V^ qc V n) )[y (n+1) - y {n) ] = r {n \ (3) 

where is the residual 



r 



W ■= / + jrqrf( 2/ (n) 



The GMRES methods proposed and analyzed in this paper apply to the solution of the 
linear Newton equations ^ or their approximations. Since the QCF equilibrium equations 
are generally solved along a quasi-static process (7], a good initial guess is usually available 
and a small number of iterations of the outer iteration ^ is sufficient to maintain stability 
and accuracy. 

1.2. Outline. We begin in Sections [2] and [3] by introducing the most important quasicon- 
tinuum approximations and outlining their stability properties, which are mostly straight- 



forward generalizations of results from [10 [12 13]. We also present careful numerical studies 
of the spectral properties of L^ cf which are particularly useful for the analysis of Krylov 
subspace methods in Section [5j 

In Section]!] we revisit the ghost force correction (GFC) scheme [27 1 which, as was pointed 
out in [6] , can be understood as a linear stationary iterative method ^ for solving the QCF 
equilibrium equations. We show that, even though the QCF method itself is stable up to a 
critical strain F*, the GFC scheme becomes unstable at a significantly reduced strain for our 
model problem. This leads us to conclude (though the simple examples we analyze here can 
only be first indicators) that the GFC method is not universally reliable near instabilities. 
We note, however, that the GFC method can be expected on the basis of both theoretical [TO] 
and computational results |10l[2l] to be more accurate near instabilities than the use of the 



uncorrected QCE energy £ qce (y), as explained in Section 4.1 Numerical results have also 



shown that the GFC method can give an accurate approximation of critical loads if the 



atomistic-to-continuum interface is sufficiently far from the defect 21, Figure 16], at a cost 



of a larger atomistic region than likely required by the accuracy of the QCF approximation. 



The quasi-nonlocal energy £ qnl (y) of (29] given by (20) is a more reliable and accurate 
energy to use in the splitting iteration d2|. It has been shown to reproduce the atomistic 
stability of one- dimensional atomistic systems with next-nearest neighbor interactions JlO] , 
and the error for multi-dimensional atomistic systems is likely to be acceptable if the longer- 
range interactions decay sufficiently fast. The splitting iteration (J2| can then be used as 
part of a continuation algorithm for a quasi-static process [7] that provides the reliable 
detection of the stability of the atomistic system [To] as well as the improved accuracy for 
the deformation given by the force-based approximation [13] . 

We conclude Section [4] by proving the numerical instability of the modified conjugate 



algorithm 21 for our QCF model problem. We present these two examples to demonstrate 
the subtleties in designing an iterative algorithm for the solution of the QCF system and 
to underscore the need for thorough numerical analysis in the development of stable and 
efficient iterative methods for the QCF system. 
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We conclude by considering in Section [5] the generalized minimal residual method (GM- 
RES) for the solution of the indefinite and non-symmetric QCF system. We provide an 
analysis of basic as well as preconditioned GMRES methods. We find in this section that 
a non-standard preconditioned GMRES method, based on the discrete jy 1,2 -inner product, 
appears to have excellent stability properties up to the critical strain and a more reliable 
termination criterion. 



2. QUASICONTINUUM APPROXIMATIONS AND THEIR STABILITY 

In this section, we give a condensed description of the prototype QC approximations and 
their stability properties. We refer the reader to 10, 12 for more details. Many details of 
this section can be skipped on a first reading and only referred back to when required. 

2.1. Notation. Before we introduce the atomistic model and its QC approximations, we 
define the notation that will be used throughout the paper. 

We consider a one-dimensional atomistic chain whose 2N + 1 atoms have the reference 
positions Xj = je for e = 1/N. We will constrain the displacement of boundary atoms which 
gives rise to the displacement space 

U = {u G R 2N+1 : U- N = u N = 0}. 

We will equip the space U with various norms which are discrete variants of the usual Sobolev 
norms that arise naturally in the analysis of elliptic PDEs. For displacements v G U and 
1 < p < oo, we define the £ p norms, 



max^=_Ar+i,...jv \ve\, p = oo, 



and we let U 0,p denote the space U equipped with the £ p norm. The inner product associated 
with the £ 2 norm is 

N 

(v, w) :— e ViWi for v, w G U. 

e=-N+i 

In fact, we use \\f\\ep and (f,g) to denote the £f-norm and £^-inner product for arbitrary 
vectors /, g which need not belong to hi. In particular, we further define the U 1,p norm 

IMIw 1 ^ := \W\\e e i ( 4 ) 

where {y')i = v' e = e~ l {vi — t = —N + 1, . . . , N, and we let U 1,p denote the space U 

equipped with the U 1,p norm. Similarly, we define the space U 2 ' p and its associated U 2,p norm, 
based on the centered second difference v" t = E~ 2 (ve +1 — 2v£ + v e _i) for £ = —N + l, . . . , iV— 1. 
(We remark that, for v G U, we have that v' G R 2N and v" G M 2iV_1 .) 

For a linear mapping A : U\ — > U2 where Ui are vector spaces equipped with the norms 
|| ■ ||^, we denote the operator norm of A 

11/111 \\ Av \\u 2 
||^||l(Wi, u 2 ) '■= SU P -n — n — • 

If Ui =U 2 , then we use the more concise notation 

||^4||wi := ||^4|U(Wi, Ui)- 
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If A : IA 0,2 — > IA 0,2 is invertible, then we can define the condition number by 

cond(A) = ||A|| w o,2 • || A~ x \\uo.2. 
When A is symmetric and positive definite, we have that 

cond(A) = \£ N _J\f 

where the eigenvalues of A are 

0<Xf <---<Xijf-v 

If a linear mapping A : IA — > IA is symmetric and positive definite, then we can define the 
v4-inner product and A-norm by 

(v, w) A ■= (Av, w), \\v\\ 2 A = (Av, v). 

We define the discrete Laplacian L : IA — )■ U by 



(Lv)j := -v) 



-v j+ i + 2vj - Vj-i 
e 2 



j = -N + l,...,N-l. (5) 



A definition of the Li 1,2 inner product and norm that is equivalent to Q can now be given 
by 

(v,w)ui,2 := (Lv,w), |Mlw*> 2 = (Lv,v) = ||L 1//2 i;||| = \\v'\\ 2 2. (6) 

Since L^ 1 : IA — > IA is symmetric and positive definite, we can also define the IA~ 1,2 inner 
product and "negative" norm by 

(v,w) u -i,2 := (L-\w), \\v\\ 2 u - ia = (L~\v) = \\L- l l 2 v\\\. (7) 

2.2. The atomistic model. We consider a one-dimensional atomistic chain whose 2N + 3 
atoms have the reference positions Xj = je for e — 1/N, and interact only with their nearest 
and next-nearest neighbors. (For an explanation why we require 2N + 3 instead of 2N + 1 
atoms as previously stated, we note that the atoms with indices ±(N + 1) will later be 
removed from the model, and refer to Remark[T]for further details.) We denote the deformed 
positions by i/j, j = —N — 1, . . . , N + 1; and we constrain the boundary atoms and their 
next-nearest neighbors to match the uniformly deformed state, = Fje, where F > is a 
macroscopic strain, that is, 

y-jv-i = ~F(N + l)e, y. N = -FNe, 

Un — FNe, y N+1 = F(N + l)e. U 

The total energy of a deformation y e IR 2Ar+3 is given by 

N 



where 



j=-N 



N+l N+l 

j=-N ' j=-N+l 

N+l N+l 

= E + E *m + vu, 



(9) 
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Here, is a scaled two-body interatomic potential (for example, the normalized Lennard- 
Jones potential, <f>(r) = r~ 12 — 2r~ 6 ), and fj, j = —N, . . . , N, are external forces. We do 
not apply a force at the atoms ±(JV + 1), which will later be removed from the model. The 
equilibrium equations are given by the force balance conditions at the unconstrained atoms, 



-^{y & ) = f 3 for j = -N + l,...,N-l, 

y) = Fje for j — —N — 1, —N, N,N+1, 

where the atomistic force (per lattice spacing e) is given by 

ld£ a (y) 



(10) 



*T(y) ■= - 



e dyj 

[4>'(y'j+i) + (p'(y' j+ 2 + y'j+i)] - [<P\y'j) + 4>'(y'j + y'j-i)] 

U±(N+1) 



)2Af+3 



2.2.1. Linearization of T & . To linearize (11) we let u G w r ", u±n 
displacement from the homogeneous state yj = Fje; that is, we define 

Vj -yf foij = -N-l,...,N + l. 



'ID 



0, be a 



U; 



We then linearize the atomistic equilibrium equations (10) about the homogeneous state y F 
and obtain a linear system for the displacement u a , 



fj 




for j = -JV + 1,...,JV- 1, 

for j = —N — 1, -N, N,N + 1, 



where (Lpv)j is given by 
{L%v) s := 0! 



+ 



'2F 



-V j+2 + 2Vj - Vj _ 2 



Here and throughout we define 



f{F) and ^ F :=<f>"{2F) t 

where <p is the interatomic potential in ([9]). We will always assume that <p" F > and fi^p < 0, 
which holds for typical pair potentials such as the Lennard- Jones potential under physically 
realistic strains F . For example, if <fi is the Lennard- Jones potential, and if 4>"(r t ) = then 
<p'(rt/2)/(j)(rt) ~ 1.2 x 10 4 . This shows that the force to compress a chain to achieve a strain 
F for which (f)"(2F) < is several orders of magnitude larger than the force to fracture the 
chain. 

2.2.2. Stability of L F . The stability properties of L F can be best understood by using a 



representation derived in 10 



N 



N 



(L F u, u) = eAp Kl 2 " £ V 2 'f \ u "\ 2 = A f\W\\% - e 2 ct)'{ F \\u"\ 



2 

a-* 



121 



e=-N+i 

where Ap is the continuum elastic modulus 

A F = c, 



-N 



" 1 A J." 
>F + 



2F- 
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Following the argument in 10, Prop. 1], we prove the following equality in 11 which 



describes the stability of the uniformly stretched chain. 
Proposition 1. If (j)'^ F < ; then 

(L F u,u) 2 „ 

mm — — -7. — = A F — e is £ (p 9 F, 

ue *™+3\{0} \\U% F ^ 2F) 

«±jv=m±(]v+i)=0 

where < v e < C for some universal constant C . 

2.2.3. The critical strain. The previous result shows, in particular, that L F is positive defi- 
nite, uniformly as iV — > oo, if and only if Ap > 0. For realistic interaction potentials, L F is 
positive definite in a ground state Fq > 0. For simplicity, we assume that Fq = 1, and we ask 
how far the system can be "stretched" by applying increasing macroscopic strains F until it 
loses its stability. In the limit as N — > oo, this happens at the critical strain F which solves 
the equation 

A Fm = 4>"(F*) + 40"(2F,) = 0. (13) 



Remark 1. We introduced the two additional atoms with indices ±(N + 1) so that the 
uniform deformation y = y F is an equilibrium of the atomistic model. As a matter of fact, 
our choice of boundary condition here is very close in spirit to the idea of "artificial boundary 



conditions" (see 13, Section 2.1] or [17]), which are normally used to approximate the effect 
of a far field. In the quasicontinuum approximations that we present next, these additional 
boundary atoms are not required. □ 

2.3. The Local QC approximation (QCL). The local quasicontinuum (QCL) approxi- 
mation uses the Cauchy-Born approximation to approximate the nonlocal atomistic model 
by a local continuum model (6,[20|24 . In our context, the Cauchy-Born approximation reads 

{e-\y l+1 - y t - x )) « \ [0(2$ + 0(2^ +1 )], 

and results in the QCL energy, for y e IR 27V+3 satisfying the boundary conditions (|8j), 

N 



-N+l 



+ e 



N 



Hv-n) + \<t>{ 2 y'-N) + <t>{y'N+\) + I^v'n+i) 



(14) 



= E eW^) + 0(2^)]+£ [20(F) + 0(2F)]. 

j=-N+l 

Imposing the artificial boundary conditions of zero displacement from the uniformly deformed 
state, = Fje, we obtain the QCL equilibrium equations 

-J^ cl ( y ^) = f j for j = -N+l,...,N-h 



yf 



Fje 



for 



-N, N, 



cS 
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where 



1 d£ qc \y) 
s dyj 



W S+ i) + 20 / (2 2/ ; +1 )] - WW + 2<p'(2y' J )} 



(15) 



In particular, we see from (15) that the QCL equilibrium equations are well-defined with 
only a single constraint at each boundary (see also Remark [I]), and we can restrict our 
consideration to y E IR 2Af+1 with the boundary conditions y^N = —F and y^ = F. 

Linearizing the QCL equilibrium equations (15) about the uniformly deformed state y F 
results in the system 

qcl 



fi 





for 
for 



-N + l, 
-N, N, 



N 



where (L F 2l v)j, for a displacement v G U, is given by 



(L 



+ 40 



" 1 

2FJ 



-v j+ i + 2vj - Vj. 



-A F v'! : 



-N + l, 



N — 1. 



The increased efficiency of the local QC approximation is obtained when its equilibrium 



equations (15) are coarsened by reducing the degrees of freedom, using piecewise linear 
interpolation between a subset of the atoms [6 20 . For the sake of simplicity of exposition, 
we do not treat coarsening in this paper. 
We note that 



Lf = A F L 



where L : U — > U is the discrete Laplacian rt5J). Since the QCL operator is simply a scaled 
discrete Laplace operator, its stability analysis is straightforward: 

(L^u, u) 



Af\\u 



l\\2 



for all u EU. 



In particular, it follows that Lp 21 is stable if and only if Ap > 0, that is, if and only if F < F*, 
where F* is the critical strain defined in (13). 



2.4. The force-based QC approximation (QCF). In order to combine the accuracy of 
the atomistic model with the efficiency of the QCL approximation, the force-based quasi- 
continuum (QCF) method decomposes the computational reference lattice into an atomistic 
region A and a continuum region C, and assigns forces to atoms according to the region they 
are located in. Since the local QC energy (14) approximates y'- + in ^ by 2y'j, it is 
clear that the atomistic model should be retained wherever the strains are varying rapidly. 
The QCF operator is given by [6l[7] 



^ Cf (y) 



JJ(y) ifj'GA 
Ff{y) if j EC, 



(16) 



N 



and the QCF equilibrium equations by 
-Ff{y^) = f 3 

vf l = F J£ 

We recall that J-" qcf is a non- conservative force field and cannot be derived from an energy [6] . 



for j = -N + 1, 
for j = -N, N. 
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For simplicity, we specify the atomistic and continuum regions as follows. We fix K £ N, 
1< K < N - 2, and define 



A = {-K, . . . , K} and C = {-N + 1, . . . , N - 1} \ A. 



Linearization of (16) about y F reads 



{Lfu^) j = f j for j = - N + l,...,N-l, 
uf = for j = -N, N, 

where the linearized force-based operator is given explicitly by 



(17) 



frqtf x ._ f (L q /v)j, for j e C, 
jj ' \ (L^'u)j, for j e A. 

We note that, since atoms near the artificial boundary belong to C, only one boundary 
condition is required at each end. 

We know from [12] that the stability analysis of the QCF operator L^ cf is highly non-trivial. 
We will therefore treat it separately and postpone it to Section [3j 



2.5. The original energy-based QC approximation (QCE). In the original energy- 
based quasicontinuum (QCE) method [24], an energy functional is defined by assigning 
atomistic energy contributions in the atomistic region and continuum energy contributions 
in the continuum region. In the context of our model problem, it can be written as 



£ qc \y) = e J2 £ tiv) + £ E E ^y) for y e 



D 2AT+1 

<-i \y) 1U1 y cu< , 



where 



£t{y) = W( 2 y'e) + M) + M+i) + Wt+i)), and 
£?(y) = + y'e) + M) + M+i) + M+i + ^+ 2 ))- 



The QCE method does not satisfy the patch test [8],[9|, 22 , 27] , which be seen from the 
existence of "ghost forces" at the interface, that is, V£ qce (|r ) = g ^ 0. Consequently, the 
linearization of the QCE equilibrium equations about y F takes the form (see [9j Section 2.4] 
and (8j Section 2.4] for more detail) 

(Lr« qce )i - 9i = fj for J = -iV + l,...,iV-l, 

uf c = for j = -N, N, [ } 
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where, for < j < N — 1, we have 



+ 



As 2 



'2F 



4- 

4- 
4- 
4- 



Ae 2 



+ 



As 2 



l Vj . 


f2 — 




e 

2 Vj . 


2e 
f i — 


v .i 


e 

2 Vj 


e 


-i 


e 


e 
~ v i 


-2 


e 


2s 





+ 



+ 



£ 2e 



•Jj-2 



2 s 



e 2 

-v j+ i + 2vj - Vj-i 

and where the vector of "ghost forces," g, is defined by 

0, 0<j<K-2, 

-Te^Fi J=K-1, 

Ye^Fi 3 = K , 

±<P' 2F , j = K + l, 

-Te<t>2F, J=K + 2, 



< j < K - 2, 
j = K-l, 
J = K, 

, J = K + 1, 
j = K + 2, 



K + 3 <j < N -1, 



9i 



0. 



K + 3 < j < N -1. 



For space reasons, we only list the entries for < j < N — 1. The equations for j = 
—N + 1, . . . , — 1 follow from symmetry. 

We prove in |TT] the following new sharp stability estimate for the QCE operator L^ ce 
which implies that the L^ ce operator gives an 0(1) approximation for the critical strain, F*. 

Lemma 2. If K > I, N > K + 2, and $ F < 0, then 

inf (L^ e u,u) = A F + Xk4>2Fi 

where | < Xk < 1- Asymptotically, as K — >■ oo, we have 

X K ~ A, + 0(e~ cK ) where A* « 0.6595 and c « 1.5826. 



This result will be used in Section 4.1 where we analyze the ghost-force correction iteration, 
interpreted as a linear stationary iterative method for L^ cf with preconditioner L^ ce . 

2.6. The quasi-nonlocal QC approximation (QNL). The QCF method is the simplest 
idea to circumvent the patch test failure of the QCE method. An alternative approach was 



suggested in 14 29 , which is based on a modification of the energy at the interface. In 
this model, a next-nearest neighbor interaction term (p(e~ l (ye+i — Vt-i)) is left unchanged if 
at least one of the atoms £ + !,£—! belong to the atomistic region or an interface region 
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(which is implicitly defined by (20)), and is otherwise replaced, preserving symmetry, by a 
Cauchy-Born approximation, 

{e-\ye +1 - y t _ x )) « \ [4>(2y' e ) + 4>(2y' e+1 )}. 

This idea leads to the energy functional 

N 

S**(y) := e £ </>($ + e £ 0(yj + y' e+1 ) +8^1 Wv'i) + W , (20) 

where we set <p(y'_ N ) = <P{i/n+i) = 0- The QNL approximation satisfies the patch test; that 
is, y = y F is an equilibrium of the QNL energy functional. 

The linearization of the QNL equilibrium equations about the uniform deformation y F is 

(Lfu^^fr for J = -N + 1,...,N-1, 

for j = -N, N, 



uf = 



where 



„ -vj+i + 2vj - vj. 



+ 



^2F 



Ae 2 

-V j+2 + 2Vj - Vj-2 _ -Vj+2 + 2v j+1 - Vj 

As 2 e 2 
^vj+i +2vj -vj-i -vj + 2vj-x - vj-2 

£ 2 £ 2 

-v j+1 + 2vj - vj^ 



< j < K - 1, 

j = K, (21) 
J=K + 1, 
K + 2 < j < N -1. 

We observe from (21 ) that L F nl is not pointwise consistent at j — K and j — K + 1. 

Repeating our stability analysis for the periodic QNL operator in [10| Sec. 3.3] verbatim, 
we obtain the following result. 

Proposition 3. If K < N — 1, and 2 f < ; then 

inf (L F u, u) = A F . 

ll«'ll/2=l 



Remark 2. Since <p2 F = (A F — <p" F )/A, the linearized operators (<fi F ) 1 L F , ((f> F ) 1 Lf l , 
(<P" F )- l Lf, (ct)" F )- l L^\ and (<P" F )- l Lf l depend only on A F /(j) F , N and K. □ 

3. Stability and Spectrum of the QCF operator 
In this section, we collect various properties of the linearized QCF operator, which are, for 



the most part, variants of our results in 12 13 . We begin by stating a result for the lack of 
positive-definiteness of L^ cf , which lies at the heart of many of the difficulties one encounters 
in analyzing the QCF method. 

Theorem 4 (Lack of Positive-Definiteness of QCF, Theorem 1, [ 13 1 ) . If(fi F > and 
4>2 F 6l \ {0} then, for sufficiently large N, the operator Lf l is not positive-definite. More 
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A F = 


0.8 




0.6 




0.4 




0.2 


0.04 


N = 8 


4.83e- 


-13 


4.26e 


-13 


3.13e 


-13 


3.41e- 


-13 


1.71e-13 


32 


1.73e- 


-11 


1.27e- 


-11 


9.55e- 


-12 


9.55e- 


-12 


1.41e-ll 


128 


8.08e- 


-10 


4.00e- 


-10 


4.07e- 


-10 


4.15e- 


-10 


4.15e-10 


512 


1.06e- 


-08 


8.73e- 


-09 


1.40e- 


-08 


8.38e- 


-09 


8.73e-09 



Table 1. The difference between the spectra of Up and Lp^ 1 . The table 
displays the £°° norm of errors in the ordered vectors of eigenvalues for various 
choices of Ap with (j)" F = 1, for increasing N, K — [v^/Vj + 1- All entries are 
zero to the precision of the eigenvalue solver. 



precisely, there exist N e N and C\ > C2 > such that, for all N > N and 2 < K < N/2, 

-CyV 1/2 < inf (Lfv,v) < -C 2 iV 1/2 . 



lb'll f 2=i 



As a consequence of Theorem |4| we analyzed the stability of in alternative norms. 
Following the proof of 12, Theorem 3] verbatim (see also [l2j Remark 3]) gives the following 
sharp stability result. 



Proposition 5. If Ap > and (f)'^ < 0, then L^ cf is invertible with 



\\(Lfy 



<l/Aj 



If Ap = 0, then L^ cf is singular. 



This result shows that 



j qcf 
F 



is operator stable up to the critical strain at which the 



atomistic model loses its stability as well (cf. Section 2.2). In the remainder of this section, 
we will investigate, in numerical experiments, the spectral properties of the L^ cf operator for 
strains F such that Ap > and ifi'^p < 0. 



3.1. Spectral properties of in U ' 2 = £ 2 . The spectral properties of the operator 
are crucial for analyzing the performance of iterative methods in Hilbert spaces. The basis 
of our analysis of L^ cf in the Hilbert space U 0,2 is the remarkable observation that, even 
though L^ cf is non-normal, it is nevertheless diagonalizable and its spectrum is identical to 
that of L F nl . We first observed this in [12| Section 4.4] for the case of periodic boundary 
conditions. Repeating the same numerical experiments for Dirichlet boundary conditions, 
we obtain similar results. Table [IJ where we display the error between the spectrum of Z/^ cf 
and L F nl , gives rise to the following conjecture. 

Conjecture 6. For all N > 4, 1 < K < N — 2, and F > 0, the operator L^ cf is 
diagonalizable and its spectrum is identical to the spectrum of L F nl . 

We denote the eigenvalues of L F nl (and i^ cf ) by 

< Xf 1 ' < -Xf < ... < Ag_ r 
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The following lemma provides a lower bound for A^ nl , an upper bound for X^n-ii and conse- 
quently an upper bound for cond(L^ nl ) = X^^i/Xf 11 . Assuming the validity of Conjecture 
m this translates directly to a result on the spectrum of L^ cf . 



Lemma 7. If K < N - 1 and 4>'i F - 0; ^ en 

Xf 11 > 2 A F , Ag_i < (A F - m F ) e~ 2 = ^e' 2 , and 



Xf 11 



•j i 



cond(L^ ) = — r < ' 1 - 



e 



Proof. It follows from Proposition [3] and (22) that 

,qnl • , (L^V, V ) ■ f (Lf l V, V) (v', V') (V f , V') 

X{ = inf — — = mf — — — > A F mf > 2A F 

{v,v) vtu (v',v') (v,v) v£U {v,v) 

since the infimum of the Rayleigh quotient (v', v')/(v , v ) is attained for v e U where vt = 
sin((N - £)tt/(2N)) (3TJ Exercise 13.9] and has the value 

inf ^1 = 4iV 2 sin 2 fM > 2 . (22) 
v&a (v,v) \4NJ - K ' 

The estimate for the maximal eigenvalue follows similarly from 

nl (Lf l v, v) 

A 2N-l — SU P — / \ 



and the representation (21). □ 



For the analysis of iterative methods, particularly the GMRES method, we are also in- 



terested in the condition number of t 
Assuming the validity of Conjecture [i 



die basis of eigenvectors of L^ cf as N tends to infinity, 
we can write L^ cf = VA qci V~ 1 where A qcf is diagonal. 
In Figure [TJ we plot the condition number for increasing values of N and K, and for vari- 
ous choices of A F with <p" F = 1 (it follows from Remark [2] that V actually depends only on 
A F /(j)' F and N). Even though it is difficult to determine from this graph whether cond(V^) is 
bounded as iV — > oo, it is fairly clear that the condition number grows significantly slower 
than log (AT). We formulate this in the next conjecture. 

Conjecture 8. Let V denote the matrix of eigenvectors for the force-based QC operator 
Lf. If A F > 0, then cond(V) = o(log(iV)) as N -»■ oo. 

3.2. Spectral properties of in U 1,2 . To study the preconditioning of L^ cf by Lp 21 = 
A F L, we consider the (generalized) eigenvalue problem 

Lfv = XLv, veU, (23) 

which can, equivalently, be written as 

L^Lfv = Xv, v G U, (24) 
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Condition Number of Eigenbasis of L"^ c 

16 



8 



o 
o 



2 



1 

8 16 32 64 128 256 512 1024 

N 

Figure 1. Condition number of the matrix V plotted against N, with atom- 
istic region size K = [\/N\ +1, and for various values of Ap, with fixed (f) F = 1. 
Here, L^ cf = VA^V^ 1 is the spectral decomposition of L^ cf . 





A F = 


0.8 




0.6 




0.4 




0.2 


0.04 


N = 8 


3.33e 


-15 


1.13e- 


-14 


1.67e- 


-15 


2.14e- 


-15 


9.99e-16 


32 


1.88e- 


-13 


1.83e- 


-13 


4.62e- 


-14 


6.48e- 


-14 


3.94e-14 


128 


1.34e- 


-12 


5.13e- 


-13 


5.72e- 


-13 


3.85e- 


-13 


5.51e-13 


512 


2.22e 


-11 


9.78e- 


-12 


7.02e- 


-12 


4.32e- 


-12 


4.56e-12 



Table 2. The difference between the spectra of L _1 L^ C and L _1 L^ n . The 
table displays the £°° norm of errors in the ordered vectors of eigenvalues for 
various choices of F, for increasing N, K = [V^J + 1, and with fixed 4> F = 1. 
All entries are zero to the precision of the eigenvalue solver. 

or as 

L-^LfL-^w = Xw, we U, (25) 

with the basis transform w = L l l 2 v, in either case reducing it to a standard eigenvalue 
problem in 

In Table [2j we display the numerical experiment that corresponds to the same experiment 
shown in Table [lj We observe that also the W 1,2 -spectra of the L^ cf and L F nl operators are 
identical to numerical precision. 

Conjecture 9. For all N > 4, 1 < K < N - 2, and F > 0, the operator L^Lf 1 is 
diagonalizable and its spectrum is identical to the spectrum of L~ 1 L F al . 
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In the following lemma we completely characterize the spectrum of L l L F al , and thereby, 
subject to the validity of Conjecture]^} also the spectrum of L~ 1 Lf i . We denote the spectrum 



ofL- 1 ^ 1 by {/if :j = l,...,2N-l}. 

Lemma 10. Let K < N — 2 and A F > 0, then the (unordered) spectrum of L^ 1 Lff L (that 
is, the U 1 ' 2 -spectrum) is given by 



qnl _ J A F - m F sin 2 , j = 1, . . . , 2K + 1, 



A i 



j = 2K + 2, ...,2N- 1. 



In particular, if tft^p < 0, then 



qnl 

maXj 

qnl 

mm,- 



1 _ 4< />2F gin 2 ( {2K+1)k 



A t 



\ 4K+4 



A F + A F 



2F ■ 2 

sin 



(4K+4) 



.4, 



+ 0{K 



-2^ 



Proof. We will use the variational representation of L^ nl from 10, Section 3.3], which reads 

K 

(Lf l u, v) = A F (u', v') - <\>' 2F e K+i ~ u e)Wt+i ~ v e) for u,veU. 

l=-K 

Summation by parts in the second term yields 

(L q F l u,v) = A F (u',v') - <j% F (Mu',v f ) for u.veU, 
where M is the 2N x 2N matrix given by 





(° 









\ 


M = 




1 -1 
-1 2 -1 

-1 2 -1 
-1 1 



















\ 








• 0/ 



and where the first and last non-zero rows are, respectively, the rows —K and K + 1. We 
call the restriction of the conjugate operator L^ nl = A F I — 4>2 F M : M. 2N — > M. 2N to the 
2N — 1 dimensional invariant gradient space Rl N = {<p e R 2N :J2i ( Pi = °) the restricted 
conjugate QNL operator L^ nl — A F I — ^ F M : IR 27V — > M 2Ar , and we note that we can write 
the eigenvalue relation (23) in weak form as 



(Lf l u, v) = (Lf'u', v') = X(u', v') Vu G U. 



(26) 



r qnl 



We can see from (26) that the 2N — 1 generalized W 1,2 -eigenvalues of Llf 1 and the standard 



o2Af 



2N. 



n2N 



are the same. If Uj are the 2N — 1 eigenvalues of L F nl 
then, letting u^' G U be the (unique) functions for which 



t -eigenvalues of L F . ^ 
with eigenvectors <pv) in 
(u^)' = ipv\ we obtain 

(LfV j \v) = (Lf\u^y,v') = usduwyy) 



Vu G hi, 
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which is equivalent to (23). 

The operator Lf l : R 2Ar R™ has a (2N -2K - 2)-multiple eigenvalue with value A F 
and corresponding orthogonal eigenvectors </? (i) € Rl N can be taken to be the projection 
onto IR* of the canonical basis vectors corresponding to the zero-diagonal entries of M. We 
will see that the remaining 2K + 1 eigenvalues of L^ nl : M% N — > take the form 

Vj = A F - 02 F ^, 

where Uj, j — 1, . . . , 2K + 1 are the non-zero eigenvalues of the non-zero block of M, which 
we denote M. It is easy to check that the eigenvectors of the matrix M are given by 

g f = cos (jir(£ + K- l/2)/(2K + 2)), t = —K, ...,K + 1, 

for j = 0, . . . , 2K + 1, and the corresponding eigenvalues by 



4 sin 2 (jn/(AK + A)), j = 0, . . . , 2K + 1. 



The first eigenvector g^°> is constant, and hence all other eigenvectors have mean zero. 
This implies that the eigenvalues Uj, j = 1, . . . , 2K + 1, give the remaining eigenvalues of 
L^ nl : IR 2iV — > M.™ . This concludes the proof of the lemma. □ 



Remark 3. Even though Lemma 



10 



gives uniform bounds on the spectrum of L^ nl in Li 1 ' 2 , 
it does not give the desired sharper result that eigenvalues are clustered, for example, at A F . 
As a matter of fact, Lemma 10 shows that this is never the case. However, we see that, if K 
remains bounded as iV — > oo, then all but a finite number of eigenvalues of L _1 / 2 L^ cf L -1 / 2 
are identically equal to A F . □ 

We conclude this study by considering the condition number of the matrix of eigenvectors 



for the eigenvalue problems (24) and (25). We write L 1 L^ cf = VA qci V 1 , where A qcf is the 



diagonal matrix of eigenvalues of L 1 Ll nl and V is the associated matrix of eigenvectors. 



In Figure |2| we have plotted numerical results for the condition number of the matrix V. 
We note that great care must be taken when computing the basis of eigenvectors since one 



eigenvalue has a high multiplicity (cf. Lemma 10). As described in Appendix |AJ the block 
structure of the matrix L~ 1 L F 2f allows us to analytically compute most of the eigenvectors 
corresponding to the high multiplicity eigenvalue and to separately compute all remaining 
eigenvectors. 

The numerical experiment displayed in Figure [2] leads to the following conjecture. 



Conjecture 11. Let V denote the matrix of eigenvectors for the preconditioned force-based 



QC operator L^Lf 1 . If A F > 0, then cond(y) = O (iV 3 ) as N -»■ oo. 

It follows from Q and Q that we can write L _1 / 2 L^ cf L -1 / 2 = WA^W' 1 where W = 
L 1 / 2 ^/ is the associated matrix of eigenvectors. In Figure |3j we have plotted numerical results 
for the condition number of the matrix W. These calculations can be simplified by observing 
that, if we define the operator D : IR 2Af_1 -»■ R 2N by Dv := v' then W T W = V T LV = 
V T D T DV. Since the condition number of a matrix A depends only on the eigenvalues of 
A T A, it follows that cond(W) = cond(W). 
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Condition Number of Eieenbasis of L 1 L C L C 

in •* 

10 



f io 6 

o 

o 

io 4 



" io 1 io 2 io 3 

N 

Figure 2. Condition number of the matrix V plotted against the system size 
N for Af/4>" f = 0.4, and various atomistic region sizes K, where L~ 1 L ( ^ i = 
VA^V' 1 is the spectral decomposition of L~ x Up . Since (</>^) _1 L^ cf depends 
only on Ap / <fi" F and N, the matrix V depends only on Ap/ <f>" F and N. For each 
curve we have cond(V^) is 0(iV 3 ), but in fact the curves appear to grow like 

jy3/2_^3/2_ 

The numerical experiment displayed in Figure [3] leads to the following conjecture. 

Conjecture 12. Let W denote the matrix of eigenvectors for the preconditioned force- based 
QC operator L'^Lf* L' 1 / 2 . If A F > 0, then cond(VT) = O (N 3 ) as N -» oo. 

4. Iterative Methods for the Nonlinear QCF System 

In this section, we briefly review and analyze two common solution methods for the QCF 
equilibrium equations. The first method, the ghost force correction (GFC) scheme, is often 
considered an independent approximation scheme rather than an iterative method for the 
solution of the QCF system. However, it was shown in |6] that the ghost force correction, 
when iterated to self-consistency, does in fact give rise to the QCF method. In the following 
section, we will show that a linearization of the GFC method predicts a lattice instability at 
a strain significantly less than the critical strain of the atomistic model. 

The second method that we discuss solves the QCF equilibrium equations by computing 
the location along the search direction where the residual is orthogonal to the search direc- 

that the indefiniteness of L^ cf implies that this method 
cannot be expected to be numerically stable for the QCF system. 

4.1. The Ghost Force Correction. After discovering that the original energy-based QC 
method (QCE) is inconsistent at the interface, a dead load correction was proposed to 
remove the so-called ghost forces [27|. The idea of this ghost force correction (GFC) is 




tion 12 11. We show in Section 4.2 
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Figure 3. Condition number of the matrix W plotted against the sys- 
tem size N for Ap/<f>" F = 0.4, and various atomistic region sizes K, where 
L-WUpL- 1 ! 2 = WA^W' 1 is the spectral decomposition of L -1 / 2 L^ cf L -1 / 2 . 
For each curve, cond(IF) is 0(N 3 ). 



the following: Since the Cauchy-Born continuum model is consistent with the atomistic 
model, the "defective" (inconsistent) forces of the QCE method at the interface are simply 
replaced by the Cauchy-Born forces in the continuum region and by the atomistic forces in 
the atomistic region. The discrepancy between the forces of the QCE method and those of 
the QCF method are called the ghost forces, and are defined as follows: 



where 



g(y):=^ c{ (y)-^ ce (y) 



It is clear that the ghost forces are concentrated in a neighborhood of the atomistic-to- 
continuum interface and can therefore be computed efficiently 27 . The GFC is then nor- 



mally applied during a quasistatic loading process. In the following example algorithm, the 
loading parameter is the macroscopic strain F > and the corresponding space of admissible 
deformations is y F = y F + U- 



GFC Iteration: 

0. Input: y<® G y x such that F**(y<®) + f 

1. For n = 1, 2, 3, . . . do 

2. Evaluate = g(y' n ~ 1 *), where y^ n ~ 

3. Find y^ G argmin {£ qce (y) — (f,y) 



=s 0; increment 5F > 

i) = y (»-i) + X 5F 
-(g (n \y):yey 1+n s F }. 
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Remark 4- Increased efficiency can be obtained by allowing nonuniform steps and multiple 
GFC iterations at a fixed load [7], thus introducing a second inner loop. For the purpose of 
the present paper, we will focus on the simpler algorithm above. □ 

We now consider the GFC iteration above for purely tensile loading which is given by 
/ = 0. We also take the initial iterate to be the uniform deformation for F — 1, that is, 
y(°) = y 1 = x, and 5F to be small. Then it is easy to see that the GFC iteration gives the 
uniform deformation y^> = y 1+n6F until 1 + nSF > F gic , where F gfc is the uniform strain at 
which L^ ce becomes unstable. We recall from Lemma [2] that L^ ce becomes unstable at F gfc 
satisfying 

where § < X K < 1 and (j)^ Fglc < 0, so F§ fc < F*. 

The critical strain F qce for the uncorrected energy £ qce (y) was investigated in |10] by 
linearizing S qce (y) about 

2/qce e argmin {£ qcc {y) : y e y F ] 

rather than about y F . It was shown, in agreement with the computational experiments in [10] 
and [21] , that the GFC method does improve the accuracy of the computation for the critical 
strain, that is, 

F qce < ^gfc < F 



See 10 for a more precise statement of these results. 



4.2. A modified conjugate gradient method. Another popular approach to solving the 
QCF equilibrium equations is to replace the univariate optimization used for step size se- 



lection in the nonlinear conjugate gradient method 23 with the computation of a step size 
where the residual is orthogonal to the current search direction [21] . More specifically, if 
is the current search direction, then this method computes y^ n+ ^~= y^> -\-a^ n 'd^ n ' such that 

(J rqci (y {n+1) ) + f, d {n) ) «0. (27) 
We can easily see that this method is numerically unstable by considering a linearization 



of (27) about the uniform configuration y F to obtain 

< - Lf + a (">tZ(»)) + /, = 0, 

or equivalently, 

-a {n) (Lf rf (n) , d {n) ) + {Lfu {n \ d^) + (/, dF>) = 0. 

However, according to Theorem [ZJ L^ cf is indefinite, which implies that there exist direc- 
tions d such that (L^d, d) = 0. Hence, if such a singular direction d is chosen (for example, 
if the initial iterate satisfies L^u^ = d) then the step size aS n ^ is undefined. More generally, 
if a direction d^ is "near" such a singular direction (for example, L'p ~ d), then the 
computation of a (n) is numerically unstable. 
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5. GMRES Solution of the Linear QCF Equations 

We now consider the generalized minimal residual method (GMRES) to find (approximate) 
solutions to the linear, force-based QC equilibrium equations 

Lfu<** = f. (28) 

GMRES is an attractive iterative method for the solution of nonsymmetric linear equations 
since the iterates satisfy a minimality property for the residual. This minimality property is 
the basis for our analysis of the convergence of the GMRES method for the solution of the 
QCF equations. 



5.1. Standard GMRES. We recall that GMRES [26] builds a sequence of Krylov subspaces 
IC m := span {r<°>, LfV°\ (L?)V°\ . . . , (Lf )™-V°>}, 

where := f — Up is the initial residual, and it finds an approximate solution 

u (m) := argmin 1 , 6u( o )+/Cm ||/ - Lf l v\\ 2 ( 29 ) 

that minimizes the ^-norm of the residual := / — Up for (28). The residual 
satisfies the minimality property 

||r( m )||^= min \\f - Lfv\\ p = min \\p m (Lf)r^\\ p (30) 

Pm(0)=l 

where 

V m = {polynomials p of degree < m}. 

It follows from (30) that depends only on n ', Ap/(fi F , N, and K. 

GMRES solves the minimization problem (29) by reducing it to a least squares problem 
for the coefficients of an l 2 £ — orthonormal sequence {t>i, . . . ,t> m +i} computed by the Arnoldi 



process. For details, see 26 , 32 



The convergence analysis does not require a symmetric matrix, and we will see that Con- 
jectures [6] and [8] regarding the spectrum of eigenvalues and conditioning of eigenvectors are 
exactly what is needed for an error analysis of GMRES applied to L^ f . 

Proposition 13. If Conjecture [6| holds, then 

'l 



1 /2A F 



< 2cond(V) ( I lk (0) lk ( 31 ) 

1 / 2Ap 



Remark 5. We recall from Conjecture M that cond(V^) = o(log(iV)). We note that the 



estimate (31) gives a reduction of the convergence rate for strains near the critical strain 
A F , = 0. 

Proof. By Conjecture [6j is diagonalizable, and we have that L^ cf = VA^V^ 1 where 
V contains the eigenvectors of L^ cf as its columns and where A qcf is the diagonal matrix of 
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eigenvalues of Lf l . We denote the set of eigenvalues of L^ f by a(Lf { ). We then have by (30) 
that 



||r (m) ||^= min \\p m {Lf)r {0) \\ p = min || V> m (A qcf )^ V 0) || 

Pm(0) = l Pm(0) = l 

< cond(y) inf llpjl f rq=fJk (0) || « 

Pm(0)=l 

Ibmll^f) = SUp |p m (A)|. 
Xea(Lf) 

By Conjecture [HJ L^ cf and L^ nl share the same spectrum, so we have that 



where 



inf ||p m || /r qcf, = inf ||p m || , T ^u < inf max |p m (A)|. 

P mev m UFmU °(L F ) PmeVm ^ m "«(L F ) - pmgPm Aqnl<A<Aqnl ™ i\ 

Pm(0) = l Pm(0)=l Pm (0)=1 1 2Ar_1 



We now recall 26 that 



inf max lp m (A)| < 2 _ 

Pm£Pm A qn '<A<A qnl V 1 + \/7 

Pm( o)=i Al £A ^-i V V I 

where 7 = l/cond(L qnl ) = X^/^Tn-v We have h Y Lemma Q that 7 < (2A F e 2 )/ '<f>" F . It thus 
follows that 



|r (m) ||^ < 2cond(V) [ \ + ^ ) ll r(0)| 



1 2 At 

<2cond(V)( ^L] || r (o)|| n 

1 + e , /M£ 



In Figures [4] and |5j we display the residual and error of the standard GMRES iterates 
when the algorithm is applied to the solution of the QCF system with right-hand side 

{1 x > 0, 
~~ ' (32) 
— 1, x < 0, 

which is smooth in the continuum region but has a discontinuity in the atomistic region. We 
also set Ap = 0.5 and <p" F — 1. We observe the slow convergence predicted by the theory of 
this section. However, we also observe alternation of slow and fast regimes, which our theory 
was unable to predict. 

5.2. Preconditioned GMRES with P = L. We next consider the GMRES algorithm left- 
preconditioned by P = L, which is the GMRES algorithm applied to the left-preconditioned 



QCF equilibrium equations 26 



L^Lfu^ = L~ l f. (33) 
We now denote the mth left-preconditioned Krylov subspace by 

t m =: span{L-V°), (L^Lf )L'V°\ . . . , (iT 1 Lf)™' 1 X"V )} 
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Application of standard GMRES to the QCF system (28) with 

1. We plot the l^-norm of the 



Figure 4. 

right-hand side (j32h , Ap = 0.5, and 



residual against the iteration number m for various choices of N and K. We 
observe the slow convergence of the residual partially predicted by the theory 
in section 5.1 We recall that there are 2N — 1 degrees of freedom. 



GMRES Error 



10 




Application of standard GMRES to the QCF system (28) with 

= 1. We plot the ^-norm of the error 
e (m) _ u (m) _ ^qcf^gaingt the iteration number m for various choices of N and 
K. We observe that ||e^ m ^||^2 closely mirrors the norm of the residual 



Figure 5. 
right-hand side (32), Ap = 0.5, and 
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and compute the minimizer 



u 



(m) ._ 



argmin„ eu(0 ) + £ m ||-L l (f - Lf v)\\ e ,. 
Proposition 14. If Conjecture [P] holds, then 

|L-V m) |L < 2cond(V) I 1 1 Hi- 1 /- 11 '" 



Vp 



er 



(34) 



Remark 6. We recall that Conjecture 11 states that cond(V) = O {N 3 ) 



Proof. As in the proof of Proposition 13 above, the residual satisfies 



| L -i r M| 



mm 
min 

Pm(0)=l 

min 

Pm(0) = l 



Pi 



L-\f-Lfv) 



(o) 



(35) 



<cond(V") inf |bm|L L -i£^ \\L V 0) | 

Pm(0) = l 



6 



10 



1 L^ cf V r is 



where V is a matrix with the eigenvectors of L _1 L^ C as its columns and V~~ X L 
the diagonal matrix A qcf . By Conjecture |9, L _1 L^ cf has the same spectrum as L _1 L^ nl , and 
by Lemma 



we have that 7 := ^T^/^n-i — ^-f/4>'f- Using the bound on the spectrum, 



we arrive at t 



A F /4>" F , and N. 



le estimate (34). It follows from (35) that L depends only on L 1 r^°\ 

□ 



Numerical experiments describing the convergence of the preconditioned GMRES method 
are displayed in Figures [6] and [7} In the first iteration, we observe a large decrease in the 
residual, which can be explained by the fact that 1 is a multiple eigenvalue. Next, we 
see that the iteration for the two cases with K = 4 converges to machine precision in 10 
iterations. This is an immediate consequence of Lemma 10, which shows that L _1 L^ cf has 
exactly 2K + 2 distinct eigenvalues. Finally, we observe precisely the convergence rate for 



the residual predicted in Proposition 14 which is independent of N and K. However, we also 
notice in Figure [7] that the error is not directly related to the residual. This may be caused 
by a large condition number of the eigenbasis, and means that the residual is not necessarily 
a reliable termination criterion. Finally, we note that, even though in this experiment Ap is 
close to zero (that is, the systems is close to an instability), we still observe rapid convergence 
of the method. 

5.3. Preconditioned GMRES with P = L in the U 1,2 norm. A possible reason for 
the poor connection between residual and error in the preconditioned GMRES method is 
that we have minimized the residual in an inappropriate norm. A more natural norm than 
||£"V m )|| £i is the W 1>2 -norm @ of L~ l r^ 



\[ J - 1 r ( rn )\ 



u 1 



|£-l/2 r (m)| 



,(m)| 



u- 
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P-GMRES: Residual 
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Figure 6. Application of preconditioned GMRES to the QCF system (|28j) 
with right-hand side (32), and with Ap = 0.1 and <p" F = 1. We plot the r e - 



norm of the preconditioned residual against the iteration number m for various 
choices of N and K. We observe precisely the convergence rate 
q m with q = (1 - A F / (j)" F ) / {I + yjAF/tf'p), predicted in Proposition [l4 



This gives a clear motivation for minimizing the preconditioned residual L~ l r^ in the U 1 ' 2 - 
norm (see also |30j Sec. 13] for a more extensive discussion of this idea and interesting 
generalizations) . 

This leads to a variant of the preconditioned GMRES method where, at the mth step, we 
compute the minimizer 



u 



\L-\f 



Lfv) 



iu 1 



by computing an Arnoldi sequence {y\, ...,v m+ i} that is U 1,2 — orthonormal for the left- 
preconditioned equations (33). We then obtain, subject to the validity of Conjecture [9j that 
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P-GMRES: Error 
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Figure 7. Application of preconditioned GMRES to the QCF system (|28j) 
with right-hand side (32), and with Ap = 0.1 and <fi F = 1. We plot the ^-norm 



of the error = u^ r ^—u qct against the iteration number m for various choices 



of N and K. The expected rate is ||e 



(m) I 



q m w here g = (l- y /A F /<//j r )/(l+ 



the residuals satisfy 

II /~-M m )|l 



= min \\ Vm (L- l Lf)L- l r^\\ uia 

I'm ev m 

Pm(0) = l 

= min II^A^^L-V )!! 

Pm(0)=l 

< cond (L 1 ^ 2 V) inf IkJI lrqcf JlL" 1 ^! 

Pm(0) = l 



u 1 



< 2 cond (W) 



i+ /^e 



I r-VCO) II 



It follows from (36) that L M" 1 ) depends only on L 1 r(°\ Ap/<p" F , and N. 
We have thus proven the following convergence result. 



Proposition 15. If Conjecture^ holds, then 
||£" lr(m) |L,a < 2 cond W 



I/— M°)H 



2(> 



M. DOBSON, M. LUSKIN, AND C. ORTNER 



Modified P-GMRES: Residual 
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Figure 8. Application of the preconditioned GMRES algorithm with U 



1.2 



inner product to the QCF system (28) with right-hand side (32), and with 
Ap = 0.1 and (j)" F = 1. We plot the £7^ 2 -norm ^ of the residual against the 
iteration number m, for various choices of iV and K. We observe precisely the 
convergence behaviour predicted by Proposition [l5| namely ||r^ m )||^-i,2 ~ q m 
where q = (1 - V^V0£)/(1 + v^W^)- 



Remark 7. We recall from Conjecture 12 that cond(W) = O (N 3 ) 



We have tested this variant of the preconditioned GMRES method for the system (28) 
with right-hand side (32) and displayed the detailed convergence behavior in Figures [8] and |9j 
All our observations about the residual that we made in the previous section are still valid; 
in particular, the spectrum of L _1 L^ nl (that is, of L _1 L^ cf ) fully predicts the convergence of 
the residual. Moreover, we notice that the residual and the error are now closely related, 
that is, the residual can be taken as a reliable termination criterion for the iterative method. 
Of course, we have not presented a proof for this statement and further investigations should 
be performed to verify this. 

To conclude we remark that, even though we find the GMRES method in the ZY 1,2 -inner 
product more attractive from a theoretical point of view, we have no evidence that it is con- 
siderably more efficient in practice than the more standard preconditioned GMRES method 
presented in Section 5.2. As a matter of fact, additional numerical experiments, the details 
of which we do not display here for space reasons, show that the decay of the error in the 
W 1,2 -norm is quite similar for both methods, for a variety of choices of N, K, and /. 



Conclusion 

We began by studying the widely used ghost force correction method (GFC) [27], which 
can be understood as a linear stationary method for QCF using the QCE operator as a 
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Modified P-GMRES:Error 
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Figure 9. Application of the preconditioned GMRES algorithm with U 



1,2 



inner product to the QCF system (28) with right-hand side (32), and with 
A F = 0.1 and <f> F = 1. We plot the W T ' 2 -norm of the error = — w qcf 
against the iteration number m, for various choices of N and K. We observe 
that ||e( m )||^i,2 closely mirrors the norm of the residual ||L _1//2 r^ m - ) ||^2, that is, 
the residual provides a reliable prediction for the actual error. 



preconditioner. We showed that the GFC method becomes unstable for our model problem 
before the critical strain is reached. In practice, this means that the ghost force correction 
method would predict a reduced critical strain for the onset of defect formation or motion. 
We also showed that a popular modified nonlinear conjugate gradient method to solve the 
QCF equations 21 is numerically unstable for our model problem. 

We then proposed and studied several variants of the generalized minimal residual method 
(GMRES), which are a natural choice for the non-symmetric QCF operator. Since our ex- 
perience with stationary methods indicates that the QCL preconditioner combines efficiency 
and reliability [11], we focused exclusively on this preconditioner. Our analysis and compu- 
tational experiments have led us to propose a GMRES method, which uses the QCL method 
as a preconditioner as well as the underlying inner product. This method is reliable for our 
model problem up to the critical strain, and the residual appears to offer a more effective 
termination criterion. 

Future research will explore the extension of the algorithms and analysis in this paper 
to the multi-dimensional and nonlinear setting to develop predictive and efficient iterative 
solution methods for more general force-based hybrid and multiphysics methods |4| |l8|[2l][28] . 
Our investigations may also prove relevant for some hybrid methods that utilize overlapping 
or bridging domains [IJ see Method III]. 
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Appendix A. Eigenbasis Computation for L' 
We note that care must be taken when computing the basis of eigenvectors since the 



eigenvalue Ap has a multiplicity of (2N — 2K — 2) (cf. Lemma 10). This renders the 



problem highly ill-conditioned and naive usage of a standard eigensolver leads to unstable 

<l>i ■ - ApLej for j = 
_1 L^ cf has the block structure 



results. 
-N + h 



To circumvent this difficulty, we observe from (\17\j that L F "ej 



-K — 3 and j = K + 3, . . . , N — 1, and hence L~ 





/ A F 




x 1 


\ 


L- l Lf = 




x 2 






\ 


x 3 


Af 


Af / 



where X 2 is a (2K + 5) x (2K + 5) matrix. From this form, we see that there are 2N — 2K — 6 
standard unit vectors that are eigenvectors corresponding to the eigenvalue Ap. According 



to Lemma 10, the multiplicity of Ap is 2N — 2K — 2, so that we have accounted for all but 
four eigenvectors of the high multiplicity eigenvalue Ap. 

Next, we reduce the dimensionality of the eigenvalue problem to 

X 2 v 2 = Xv 2 . 

We then extend these eigenvectors to eigenvectors of L _1 L^ cf by defining 

v 2 

where vi := (X — Ap)~ 1 Xiv 2 and V3 := (A — Af)~ 1 X3V 2 . Note that Vi (i = 1,3) is well defined 
provided that A 7^ Ap or XiV 2 = 0, and we observe numerically that XiV 2 = whenever 
X = Ap. Finally, the eigenvectors obtained in this manner are normalized before computing 
the condition number of the eigenbasis. 
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